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Q ' Abstract 

CN ■ Both orbital and attitude dynamics employ the method of variation of parameters. In 

a non-perturbed setting, the coordinates (or the Euler angles) get expressed as functions 
of the time and six adjustable constants called elements. Under disturbance, each such 
" ' expression becomes ansatz, the "constants" being endowed with time dependence. The 
. perturbed velocity (linear or angular) consists of a partial time derivative and a convective 
I I term containing time derivatives of the "constants." It can be shown that this construc- 
^ ■ tion leaves one with a freedom to impose three arbitrary conditions upon the "constants" 
. and/or their derivatives. Out of convenience, the Lagrange constraint is often imposed. 
^ I It nullifies the convective term and thereby guarantees that under perturbation the func- 
' tional dependence of the velocity upon the time and "constants" stays the same as in the 

■ undisturbed case. "Constants" obeying this condition are called osculating elements. 

I The "constants" chosen to be canonical, are called Delaunay elements, in the orbital 

' case, or Andoyer elements, in the spin case. (As some of the Andoyer elements are time- 
p i' dependent even in the free-spin case, the role of "constants" is played by these elements' 
Q , initial values.) The Andoyer and Delaunay sets of elements share a feature not readily 
I apparent: in certain cases the standard equations render these elements non-osculating. 

■ In orbital mechanics, elements calculated via the standard planetary equations come 
• • , out non-osculating when perturbations depend on velocities. To keep elements osculating 

. I under such perturbations, the equations must be amended with extra terms that are not 
^ ' parts of the disturbing function (Efroimsky and Goldreich 2003, 2004). For the Kepler el- 
. ements, this merely complicates the equations. In the case of Delaunay parameterisation, 
these extra terms not only complicate the equations, but also destroy their canonicity. 
So under velocity-dependent disturbances, osculation and canonicity are incompatible. 

Similarly, in spin dynamics the Andoyer elements come out non-osculating under 
angular-velocity-dependent perturbation (a switch to a noninertial frame being one such 
case). Amendment of the dynamical equations only with extra terms in the Hamiltonian 
makes the equations render nonosculating Andoyer elements. To make them osculating, 
more terms must enter the equations (and the equations will no longer be canonical). 

It is often convenient to deliberately deviate from osculation by substituting the La- 
grange constraint with an arbitrary condition that gives birth to a family of nonosculat- 
ing elements. The freedom in choosing this condition is analogous to the gauge freedom. 
Calculations in nonosculating variables are mathematically valid and sometimes highly 
advantageous, but their physical interpretation is nontrivial. For example, nonosculat- 
ing orbital elements parameterise instantaneous conies not tangent to the orbit, so the 
nonosculating inclination will be different from the real inclination of the physical orbit. 

We present examples of situations in which ignoring of the gauge freedom (and of the 
unwanted loss of osculation) leads to oversights. 
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1 Introduction 



1.1 Historical Prelude 

The orbital dynamics is based on the variation-of-parameters method, invention whereof is 
attributed to Euler (1748, 1753) and Lagrange (1778, 1783, 1808, 1809, 1810). Though both 
greatly contributed to this approach, its initial sketch was offered circa 1687 by Newton in his 
unpublished Portsmouth Papers. Very succinctly, Newton brought up this issue also in Cor. 3 
and 4 of Prop. 17 in the first book of his Principia. 

Geometrically, the part and parcel of this method is representation of an orbit as a set of 
points, each of which is contributed by a member of some chosen family of curves C{k), where 
K stands for a set of constants that number a particular curve within the family. (For example, 
a set of three constants k — {a, b, c} defines one particular hyperbola y — ax^ + bx + c 
out of many). This situation is depicted on Fig.l. Point A of the orbit coincides with some 
point Ai on a curve C{k,i). Point B of the orbit coincides with point A2 on some other curve 
C{k2) of the same family, etc. This way, orbital motion from A to B becomes a superposition 
of motion along from Ai to A2 and a gradual distortion of the curve from the shape 
C{k,i) to the shape C{k2). In a loose language, the motion along the orbit consists of steps 
along an instantaneous curve C{k) which itself is evolving while those steps are being made. 




Fig.l. Each point of the orbit is contributed by a member of some family of 
curves C{k) of a certain type, k standing for a set of constants that number 

a particular curve within the family. Motion from A to B is, first, due to the 
motion along the curve C{k) from Ai to A2 and, second, due to the fact 
that during this motion the curve itself was evolving from C{ki) to C{k2) ■ 
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Normally, the family of curves is chosen to be that of ellipses or that of hyperbolae, k being 
six orbital elements, and A being the time. However, if we disembody this idea of its customary 
implementation, we shall see that it is of a far more general nature and contains three aspects: 

1. A trajectory may be assembled of points contributed by a family of curves of an essen- 
tially arbitrary type, not just conies. 

2. It is not necessary to choose the family of curves tangent to the orbit. As we shall see 
below, it is often beneficial to choose them nontangent. We shall also see examples when in 
orbital calculations this loss of tangentiality (loss of osculation) takes place and goes unnoticed. 

3. The approach is general and can be applied, for example, to Euler's angles. A disturbed 
rotation can be thought of as a series of steps (small turns) along different Eulerian cones. An 
Eulerian cone is an orbit (on the Euler angles' manifold) corresponding to an unperturbed spin 
state. Just as a transition from one instantaneous Keplerian conic to another is caused by dis- 
turbing forces, so a transition from one instantaneous Eulerian cone to another is dictated by 
external torques or other perturbations. Thus, in the attitude mechanics, the Eulerian cones 
play the same role as the Keplerian conies do in the orbital dynamics. Most importantly, a 
perturbed rotation may be "assembled" of the Eulerian cones in an osculating or in a nonoscu- 
lating manner. An unwanted loss of osculation in attitude mechanics happens in the same 
way as in the theory of orbits, but is much harder to notice. On the other hand, a deliberate 
choice of nonosculating rotational elements in attitude mechanics may sometimes be beneficial. 

From the viewpoint of calculus, the concept of variation of parameters looks as follows. We 
have a system of differential equations to solve ("system in question") and a system of differ- 
ential equations ("fiducial system") whose solution is known and contains arbitrary constants. 
We then use the known solution to the fiducial system as an ansatz for solving the system in 
question. The constants entering this ansatz are endowed with time dependence of their own, 
and the subsequent substitution of this known solution into the system in question yields equa- 
tions for the "constants." The number of "constants" often exceeds that of equations in the 
system to solve. In this case we impose, by hand, arbitrary constraints upon the "constants." 
For example, in the case of a reduced A^-body problem, we begin with 3(A^ — 1) unconstrained 
second-order equations for 3{N — 1) Cartesian coordinates. After a change of variables from 
the Cartesian coordinates to the orbital parameters, we end up with 3{N — 1) differential 
equations for the Q(N — 1) orbital variables. Evidently, 3(A^ — 1) constraints are necessary.^ 
To this end, the so-called Lagrange constraint (the condition of the instantaneous conies being 



^ In a fixed Cartesian frame, any solution to the unperturbed reduced 2-body problem can be written as 

fj{t, Ci , . . . , Ce) , J 
9 jit, Ci, . . . , Ce) , gj 



Xj = fj{t, Ci, . . . , Ce) , i = 1 , 2 , 3 , 

^tJc 

the adjustable constants C standing for orbital elements. Under disturbance, the solution is sought as 



xj = fj{t, C,{t), Ce{t)) , i = 1, 2, 3 , 

XJ = gj{t, Ci{t), Ce{t)) + ^j{t, C,{t), Ce{t)) , gj ^ f^)^ , ^ ^ c_ 



Insertion of Xj = fj{t. C) into the perturbed gravity law yields three scalar equations for six functions Cr{t). 
This necessitates imposition of three conditions upon Cr and Cr- Under the simplest choice = , j = 1,2,3, 
the perturbed physical velocity Xj{t, C) has the same functional form as the unperturbed gj{t, C) . Therefore, 
the instantaneous conies become tangent to the orbit (and the orbital elements Cr{t) are called osculating). 
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tangent to the physical orbit) is introduced almost by default, because it is regarded natural. 
Two things should be mentioned in this regard: 

First, what seems natural is not always optimal. The freedom of choice of the supplementary 
condition (the gauge freedom) gives birth to an internal symmetry (the gauge symmetry) of 
the problem. Most importantly, it can be exploited for simplifying the equations of motion for 
the "constants." On this issue we shall dwell in the current paper. 

Second, the entire scheme may, in principle, be reversed and used to solve systems of differ- 
ential equations with constraints. Suppose we have N + M variables Cj{t) obeying a system 
of differential equations of the second order and M constraints expressed with first-order 
differential equations or with algebraic expressions. One possible approach to solving this sys- 
tem will be to assume that the variables Cj come about as constants emerging in a solution to 
some fiducial system of differential equations. Then our second-order differential equations 
for Cj{t) will be interpreted as a result of substitution of such an ansatz into the fiducial 
system with some perturbation, while our M constraints will be interpreted as weeding out of 
the redundant degrees of freedom. Unfortunately, this subject is out of the scope of our paper, 
and it will be developed somewhere else. 



1.2 The simplest example of gauge freedom. 

Variation of constants first emerged in the nonlinear context of celestial mechanics and later 
became a universal tool. We begin with a simple example offered in Newman & Efroimsky 
(2003) 

A harmonic oscillator disturbed by a force AF(t) gives birth to the initial-condition problem 



X -\- X — ^F{t) , with x{0) and x{0) known , 
whose solution may be sought using ansatz 

X — Ci{t) sin t + C2{t) cost . 

This will lead us to 

X ^\Ci{t) sin t + C2{t) cost] + Ci{t) cost - C2{t) sin t 



(1) 



(2) 



(3) 



It is common, at this point, to put the sum Ci{t) sint + C2{t) cost equal to zero, in order to 
remove the ambiguity stemming from the fact that we have only one equation for two variables. 
Imposition of this constraint is convenient but not obligatory. A more general way of fixing 
the ambiguity may be expressed as 



Ci{t) sint + C2{t) cost = (j){t) , 
0(t) being an arbitrary function of time. This entails: 

X — 4> + Ci{t) cost — C2{t) sint — Ci{t) sint — C2{t) cost , 
summation whereof with (2) gives: 

X + X — (j) + Ci{t) cost — C2{t) sint . 



(4) 



(5) 



(6) 
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Substitution thereof into (1) yields the dynamical equation re-written in terms of the "con- 
stants" Ci , C2 ■ This equation, together with identity (4), will constitute the following system: 



This leads to 



+ Ci{t) cost - C2{t) sini = AF(t) , 
Ci{t) suit + C2{t) cost = (t){t) , 



Ci = AF cost - ^ {(pcost) 

(JjL 



C2 = -AF sint + (0sint) , 

dt 



(7) 



(8) 



the function (j)(t) still remaining arbitrary.^ Integration of (8) entails: 

Ci = J cost' dt' - (pcost + oi 



(9) 



C2 ^ - J AF sint' dt' + 0sint + 02 . 
Substitution of (9) into (2) leads to complete cancellation of the terms: 

X = Ci sin t + C2 cos t = — cost j AF sint' dt' + sint j AF cost' dt' + Oi sint + 02 cost (10) 

Naturally, the physical trajectory x{t) remains invariant under the choice of gauge function 
0(t) , even though the mathematical description (9) of this motion in terms of the parameters 
C is gauge dependent. It is, however, crucial that a numerical solution of the system (8) will 
come out 0- dependent, because the numerical error will be sensitive to the choice of 0(t). This 
issue is now being studied by P. Gurfil and I. Klein, and the results are to be published soon. 
(Gurfil & Klein 2006) 

It remains to notice that (8) is a simple analogue to the Lagrange-type system of planetary 
equations, system that, too, admits gauge freedom. (See subsection 2.2 below.) 



1.3 Gauge freedom under a variation of the Lagrangian 

The above example permits an evident extension. (Efroimsky 2002a,b.) Suppose some me- 
chanical system obeys the equation 

r = F(t, r, r) , (11) 

whose solution is known and has a functional form 

r ^ fi t, Ci, Ce) , (12) 

^ Function (p{t) can afford being arbitrary, no matter what the initial conditions are to be. Indeed, for fixed 
a;(0) and x{0), the system C2(0) = a;(0) , 0(0) + Ci(0) = x{0) solves for Ci(0) and C2(0) for an arbitrary 
choice of 0(0). 
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Cj being adjustable constants to vary only under disturbance. 

When a perturbation Ai*" gets switched on, the system becomes: 

r = F{t, r, r) + AF {t , r, r) , (13) 

and its solution will be sought in the form of 

r = fit, C,{t), Ce{t)) . (14) 



Evidently, 



In defiance of what the textbooks advise, we do not put $ nil. Instead, we proceed further to 

dot standing for a full time derivative. If we now insert the latter into the perturbed equation 
of motion (13) and if we recall that, according to (11),^ d^f/dt^ — F , then we shall obtain 
the equation of motion for the new variables Cj {t) : 

E TST^T- + « = (17) 

where 



dt dCj 



so far is merely an identity. It will become a constraint after we choose a particular functional 

df • 

form $ (i ; Ci , . . . , Ce) for the gauge function $ , i.e., if we choose that the sum Cj 

be equal to some arbitrarily fixed function $(t ; Ci , .... Cg) of the time and of the variable 
"constants." This arbitrariness exactly parallels the gauge invariance in electrodynamics: on 
the one hand, the choice of the functional form of $(t ; Ci , . . . , Cq) will never^ influence the 
eventual solution for the physical variable r ; on the other hand, though, a qualified choice 
may considerably simplify the process of finding the solution. To illustrate this, let us denote 
by g(t , Ci , . . . , Cq) the functional dependence of the unperturbed velocity on the time and 
adjustable constants: 

g(t, Ci, Cq) ^ ^/(t, Ci, Co) , (19) 



^ We remind that in (11) there was no difference between a partial and a full time derivative, because at 
that point the integration "constants" C, were indeed constant. Later, they acquired time dependence, and 
therefore the full time derivative implied in (15 - 16) became different from the partial one implied in (11). 

Our usage of words "arbitrary" and "never" should be limited to the situations where the chosen gauge 
(21) does not contradict the equations of motion (20). This restriction, too, parallels a similar one present in 
field theories. Below we shall encounter a situation where this restriction becomes crucial. 
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and rewrite the above system as 



E ^ Q = - * + (20) 



S^C,^*. (21) 

If we now dot-multiply the first equation with Of /dCi and the second one with dg,/dCi , 
and then take the difference of the outcomes, we shall arrive at 

i:iC„QlC, = (Ai.-*).^-*.|| , (22) 



the Lagrange brackets being defined in a gauge-invariant (i.e., ^>-independent) fashion.^ If 
we agree that $ is a function of both the time and the parameters C„ , but not of their 
derivatives,® then the right-hand side of (22) will implicitly contain the first time derivatives 
of Cn ■ It will then be reasonable to move these to the left-hand side. Hence, (22) will be 
reshaped into 

^ (f^'^^^] + ^ ■ ^ ) ^ ^ ^ ■ " ^ ■ ^ " ^ ■ * • ^^^^ 

This is the general form of the gauge-invariant perturbation equations, that follows from 
the variation-of-parameters method applied to problem (13), for an arbitrary perturbation 
F{r, r, t) and under the simplifying assumption that the arbitrary gauge function $ is 
chosen to depend on the time and the parameters C„ , but not on their derivatives.^ As- 
sume that our problem (13) is not simply mathematical but is an equation of motion for some 

^ The Lagrange-bracket matrix is defined in a gauge-invariant way: 

3 J J 



and so is its inverse, the matrix composed of the Poisson brackets 

ir r \ = ^ ^ _ ^ ^ 



Evidently, (22) yields 



^ The necessity to fix a functional form of ^{t; C\ , . . . , Ce ) , i.e., to impose three arbitrary conditions 
upon the "constants" Cj , evidently follows from the fact that, on the one hand, in the ansatz (14) we have six 
variables C„(t) and, on the other hand, the number of scalar equations of motion (i.e., Cartesian projections 
of the perturbed vector equation (13) ) is only three. This necessity will become even more mathematically 
transparent after we cast the perturbed equation (13) into the normal form of Cauchy. (See Appendix.) 

^ We may also impart the gauge function with dependence upon the parameters' time derivatives of all 
orders. This will yield higher-than-first-order derivatives in equation (23). In order to close this system, one 
will then have to impose additional initial conditions, beyond those on r and r . 
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physical setting, so that is a physical force corresponding to some undisturbed Lagrangian 
Lo , and Ai^ is a force perturbation generated by a Lagrangian variation A£ . If, for exam- 
ple, we begin with Co{r , r , t) = r ^/2 — U{r , t) , momentum p = r , and Hamiltonian 
Ho{r , p , t) = + U{r , t) , then their disturbed counterparts will read: 

C{r,r,t) ^ - U{r) + A£(r, r, t) , (24) 

P - r + -Qf- ' (25) 
n^pr-C^^ + U + An , (26) 

The Euler- Lagrange equations written for the perturbed Lagrangian (24) are: 
where the disturbing force is given by 

Its substitution in (23) yields the generic form of the equations in terms of the Lagrangian 
disturbance (Efroimsky & Goldreich 2004): 



d (dAC \ \ dC. 



dCj \ dr J J dt 



(30) 



d 



dg _ df d _ dAC d \ ^ dAC' 



dCn dCr, dt dr dC„ \ dr 



This equation not only reveals the convenience of the special gauge 

(which reduces to # = in the case of velocity-independent perturbations) , but also explicitly 
demonstrates how the Hamiltonian variation comes into play: it is easy to notice that, according 
to (27), the sum in square brackets on the right-hand side of (30) is equal to — AH , so the 
above equation takes the form Y^j [Cn Cj] Cj — — dATi/dCn . All in all, it becomes 
clear that the trivial gauge, ^ — , leads to the maximal simplification of the variation-of- 
parameters equations expressed through the disturbing force: it follows from (22) that 

df 

^[CnCj]Cj = AF ■ — ^ , provided we have chosen * = . (32) 

dCn 
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However, the choice of the special gauge (31) entails the maximal simplification of the variation- 
of-parameters equations when they are formulated via a variation of the Hamiltonian: 

[Cn Cj] = — ~Q^~' > provided we have chosen $ = — ~~q^ ■ (33) 

It remains to spell out the already obvious fact that, in case the unperturbed force F is 
given by the Newton gravity law (i.e., when the undisturbed setting is the reduced two-body 
problem) , then the variable "constants" C„ are merely the orbital elements parameterising a 
sequence of instantaneous conies out of which we "assemble" the perturbed trajectory through 
(14). When the conies' paramctcrisation is chosen to be via the Kepler or the Dclaunay vari- 
ables, then (30) yields the gauge-invariant version of the Lagrangc-type or the Delaunay-type 
planetary equations, accordingly. Similarly, (22) implements the gauge-invariant generalisation 
of the planetary equations in the Euler-Gauss form. 

From (22) we sec that the Euler-Gauss-type planetary equations will always assume their 
simplest form (32) under the gauge choice $ = . In astronomy this choice is called 
"the Lagrange constraint." It makes the orbital elements osculating, i.e., guarantees that the 
instantaneous conies, parameterised by these elements, are tangent to the perturbed orbit. 

Prom (33) one can easily notice that the Lagrange- and Delaunay-type planetary equations 
simplify maximally under the condition (31). This condition coincides with the Lagrange 
constraint $ = when the perturbation depends only upon positions (not upon velocities or 
momenta). Otherwise, condition (31) deviates from that of Lagrange, and the orbital elements 
rendered by equation (33) are no longer osculating (so that the corresponding instantaneous 
conies are no longer tangent to the physical trajectory). 

Of an even greater importance will be the following observation. If we have a velocity- 
dependent perturbing force, we can always find the appropriate Lagrangian variation and, 
therefrom, the corresponding variation of the Hamiltonian. If now we simply add the negative 
of this Hamiltonian variation to the disturbing function, then the resulting equations (33) will 
render not the osculating elements but orbital elements of a different type, ones satisfying 
the non-Lagrange constraint (31). Since the instantaneous conies, parameterised by such non- 
osculating elements, will not be tangent to the orbit, then physical interpretation of such 
elements may be nontrivial. Besides, they will return a velocity different from the physical 
one.^ This pitfall is well camouflaged and is easy to fall in. 

These and other celestial-mechanics applications of the gauge freedom will be considered 
in detail in section 2 below. 



1.4 Canonicity versus osculation 

One more relevant development will come from the theory of canonical perturbations. Suppose 
that in the absence of disturbances we start out with a system 

. _ m^"^ . _ _ dn^°^ 

^ dp ^ ^ dq 

q and p being the Cartesian or polar coordinates and their conjugated momenta, in the 
orbital case, or the Eulcr angles and their momenta, in the rotation case. Then we switch, via 



^ We mean that substitution of the values of these elements in g{t ; Ci{t) , . . . , Ce{t)) will not give the 
right velocity. The correct physical velocity will be given by r = g -|- * . 
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a canonical transformation 

q = f{Q, P,t) , p = xiQ, P, t) (35) 

to 

Q-iy* fn-/* 
Q=^=0 , P=-_ = . «-=0, (36) 

where Q and P denote the set of Delaunay elements, in the orbital case, or the initial values 
of the Andoyer variables, in the case of rigid-body rotation. 

This scheme relies on the fact that, for an unperturbed motion (i.e., for an unperturbed 
Keplerian conic, in an orbital case; or for an undisturbed Eulerian cone, in the spin case) a 
six-constant parameterisation may be chosen so that: 

1. the parameters are constants and, at the same time, are canonical variables {Q , P} 
with a zero Hamiltonian H*{Q, P) — 0; 

2. for constant Q and P , the transformation equations (35) are mathematically equiva- 
lent to the dynamical equations (34). 

Under perturbation, the "constants" Q, P begin to evolve so that, after their substitution into 
q = fiQit), Pit),t) , p = xiQit) , Pit) , t) , (37) 
(/, X being the same functions as in (35) ), the resulting motion obeys the disturbed equations 

We also want our "constants" Q and P to remain canonical and to obey 

• _ d{n* + An*) . _ d{n* + An*) 

^ - dP ' ^ - dQ ^^^^ 

where 

7i* = and An*{Q, P t) ^ An{q{Q,P,t), p{Q,P,t), t) . (40) 

Above all, we wish the perturbed "constants" C — Q, P (the Delaunay elements, in the 
orbital case; or the initial values of the Andoyer elements, in the spin case) to osculate. This 
means that we want the perturbed velocity to be expressed by the same function of Cj{t) and 
t as the unperturbed velocity. Let us check when this is possible. The perturbed velocity is 

q ^ g + ^ (41) 

where 
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is the functional expression for the unperturbed velocity, while 



HC{t), t)^± djmiA c,it) (43) 

is the convective term. Since we chose the "constants" Cj to make canonical pairs {Q, P) 
obeying (39 - 40), then insertion of (39) into (43) will result in 

h dQn ^"^^ ^ k dPn ^ dp ■ ^^^^ 

So canonicity is incompatible with osculation when A7i depends on p. Our desire to keep 
the perturbed equations (39) canonical makes the orbital elements Q , P nonosculating in 
a particular manner prescribed by (44). This breaking of gauge invariance reveals that the 
canonical description is marked with "gauge stiffness" (term suggested by Peter Goldrcich). 

We see that, under a momentum-dependent perturbation, we still can use the ansatz (37) for 
calculation of the coordinates and momenta, but can no longer use q = dq/dt for calculating 
the velocities. Instead, we must use q = dq/dt + dATi/dp, and the elements Cj will no 
longer be osculating. In the case of orbital motion (when Cj are the nonosculating Delaunay 
elements), this will mean that the instantaneous ellipses or hyperbolae parameterised by these 
elements will not be tangent to the orbit. (Efroimsky & Goldreich 2003.) In the case of spin, 
the situation will be similar, except that, instead of an instantaneous Kcplcrian conic, one 
will be dealing with an instantaneous Eulerian cone - a set of trajectories on the Euler-angles 
manifold, each of which corresponds to some non-perturbed spin state. (Efroimsky 2004.) 

The main conclusion to be derived from this example is the following: whenever we en- 
counter a disturbance that depends not only upon positions but also upon velocities or mo- 
menta, implementation of the afore described canonical-perturbation method necessarily yields 
equations that render nonosculating canonical elements. It is possible to keep the elements os- 
culating, but only at the cost of sacrificing canonicity. For example, under velocity-dependent 
orbital perturbations (like inertial forces, or atmospheric drag, or relativistic correction) the 
equations for osculating Delaunay elements will no longer be Hamiltonian (Efroimsky 2002a,b). 

Above in this subsection we discussed the disturbed velocity q . How about the disturbed 
momentum? For sufficiently simple unperturbed Hamiltonians, it can be written down very 
easily. For example, for H — Ho + AH — p^ /2m -\- U{q) + AH we get: 

dAC ^ dAC / dAH\ ,,,, 

In this case, the pcrtTirbed momentum p coincides with the unperturbed one, g. In application 
to the orbital motion, this means that contact elements (i.e., the nonosculating orbital elements 
obeying (31) ), when substituted in g{t ; Ci , . . . , Ce), furnish not the correct perturbed velocity 
but the correct perturbed momentum, i.e., they osculate the orbit in phase space. Existence 
of such elements was pointed out long ago by Goldreich (1965) and Brumberg et al (1971). 
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2 Gauge freedom in the theory of orbits. 



2.1 Geometrical meaning of the arbitrary gauge function $ 

As explained above, the content of subsection 1.3 becomes merely a formulation of the Lagrange 
theory of orbits, provided F stands for the Newton gravity force, so that the undisturbed 
setting is the two-body problem. Then (22) expresses the gauge-invariant (i.e., taken with 
an arbitrary gauge ^{t ; Ci , . . . , Cq) ) planetary equations in the Euler-Gauss form. These 
equations render orbital elements that are, generally, not osculating. Equation (32) stands for 
the customary Euler- Gauss-type system for osculating (i.e., obeying $ = 0) orbital elements. 

Similarly, equation (30) stands for the gauge-invariant Lagrange-type or Delaunay-type 
(dependent upon whether stand for the Kepler or Delaunay variables) equations. Such 
equations yield elements, which, generally, are not osculating. In those equations, one could 
fix the gauge by putting # = , thus making the resulting orbital elements osculating. How- 
ever, this would be advantageous only in the case of velocity-independent AC . Otherwise, a 
maximal simplification is achieved through a deliberate refusal from osculation: by choosing 
the gauge as (31) one ends up with simple equations (33). Thus, gauge (31) simplifies the 
planetary equations. (See equations (46 - 57) below.) Besides, in the case when the Delau- 
nay parameterisation is employed, this gauge makes the equations for the Delaunay variables 
canonical for reasons explained above in subsection 1.4. 

The geometrical meaning of the convective term # becomes evident if we recall that a 
perturbed orbit is assembled of points, each of which is donated by one representative of a 
sequence of conies, as on Fig.2 and Fig.3 where the "walk" over the instantaneous conies may 
be undertaken either in a non-osculating manner or in the osculating manner. The physical 
velocity r is always tangent to the perturbed orbit, while the unperturbed Keplcrian velocity 
g = df /dt is tangent to the instantaneous conic. Their difference is the convective term 
So if we use non-osculating orbital elements, then insertion of their values in f(t ; Ci , . . . , Cq) 
will yield a correct position of the body. However, their insertion in g{t ; Ci , . . . , Cq) will 
NOT give the right velocity. To get the correct velocity, one will have to add ^ . (See Appendix 
for a more formal mathematical treatment in the normal form of Gauchy.) 

When using non-osculating orbital elements, we must always be careful about their physical 
interpretation. On Fig. 2, the instantaneous conies are not supposed to be tangent to the 
orbit, nor are they supposed to be even coplanar thereto. (They may be even perpendicular 
to the orbit! - why not?) This means that, for example, the non-osculating element i may 
considerably differ from the real, physical inclination of the orbit. 

We would add that the arbitrariness of choice of the function ^{t , Ci{t) , . . . , CQ{t)) had 
been long known but never used in astronomy until a recent effort undertaken by several authors 
(Efroimsky 2002a,b ; Newman & Efroimsky 2003; Slabinski 2003; Efroimsky & Goldreich 
2003, 2004; Gurfil 2004; Efroimsky 2005c) Substitution of the Lagrange constraint $ = 
with alternative choices does not influence the physical motion, but alters its mathematical 
description (i.e., renders different values of the orbital parameters Ci{t)). Such invariance of a 
physical picture under a change of parameterisation goes under the name of gauge freedom. It 
is a part and parcel of electrodynamics and other field theories. In mathematics, it is described 
in terms of fiber bundles. A clever choice of gauge often simplifies solution of the equations 
of motion. On the other hand, the gauge invariance may have implications upon numerical 
procedures. We mean the so-called "gauge drift," i.e., unwanted displacement in the gauge 
function caused by accumulation of numerical errors in the constants. 



12 



Fig. 2. The orbit is a set of points, each of which is donated by one of the 
confocal instantaneous elhpses that are not supposed to be tangent or even 
coplanar to the orbit. As a result, the physical velocity r (tangent to the orbit) 
differs from the unperturbed Keplerian velocity g (tangent to the ellipse). To 
parameterise the depicted sequence of non-osculating ellipses, and to single it 
out of the other sequences, it is suitable to employ the difference between r and 
g, expressed as a function of time and six (non-osculating) orbital elements: 
*(t , Ci , . . . , Ce) = r(t , Ci , . . . , Ce) - g(t , , . . . , Cg) . Evidently, 

dr dCj ■ ^ 

^= + = ^ + * ' 

where the unperturbed Keplerian velocity is g = dr/dt. The convective term, 
which emerges under perturbation, is $ = (dr/dCj) Cj. When a particular 
functional dependence of $ on time and the elements is fixed, this function, 
*(t , Ci , . . . , Ce), is called gauge function or gauge velocity or, simply, gauge. 




Fig. 3. The orbit is represented by a sequence of confocal instantaneous ellipses 
that are tangent to the orbit, i.e., osculating. Now, the physical velocity r 
(tangent to the orbit) coincides with the unperturbed Keplerian velocity g 
(tangent to the ellipse), so that their difference $ vanishes everywhere: 

*(t, Cu ...,Ce)^ rit,Ci, ...,Ce)-g{t,Ci, ...,Ce)=Y.-g^Cj = 0. 

This equality, called Lagrange constraint or Lagrange gauge, is the necessary 
and sufficient condition of osculation. 
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2.2 Gauge-invariant planetary equations 
of the Lagrange and Delaunay types 

We present the gauge-invariant Lagrange- and Delaunay-type equations, following Efroimsky & 
Goldreich (2003). These equations follow from (30) if we take into account the gauge- invariance 
(i.e., the ^-independence) of the Lagrange-bracket matrix [CjCj] . 



da 
'di 



2 

na 



d{-An) dAC d 



* + 



dr 9M„ 



dA£\ dg 



* + 



dAC' 
dr 



dr J dMo dMo dt 



Of d ( ^ dAC: 
dr 



(46) 



de 
di 



na? e 



d(-An) dAC d ^ dAC 



dr da 



dr 



'^ + ^]^ ^^^^ 



dr J dMo dMo dt 



dr 



(47) 



(1 - e') 
na^ e 



2U/2 



d{-An) dAC d 



duj 



dr duj 



dAC: 
dr 



^^^dAC\ds_dll/^dAC\ 



dr J duj duj dt 



dr ) 



dou — cos i 

dt na^(l — e^y/"^ sin z 



d{-An) dACd_(^dAC' 



di 



dr di 



dr 



^^^d_AC\d_^_d_l^(^d_AC: 
, dr I di di dt \ dr 



+ 



(48) 



(l-e^)V^ 
na^ e 



d(-An) dAC d ^ dAC 



de 



dr de 



dr 



^^^dAC^d^ _df d^(^^ dAC^ 



dr I de de dt 



dr 
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di 
di 



cos I 



na^ (1 — e^)-*^/^ sin i 



dj-AH) 



dAC d / dAC 

7^1^ + 



dr duj 



df 



* + 



dAC\ dg 
dr J duj 



dld_(^dAC\ 
duj dt\ dr J 



na^ (1 — e^)^/^ sin i 



d{-An) 9AC_d_f^ 



dr dil 



dAC' 
dr 



(49) 



'^^dAC\dg _dll( dAC\- 
^ dr J dn dn dt\ dr ) 



dn 



no? {1 — e^yl"^ si 



sm I 



d{-AH) dAC d , dAC 



di 



dr di 



dr 



'^^dAC\d_^_d_ld^f^d_AC: 
^ dr j di di dt\ dr ^ 



(50) 



dMo 
"df 



na^ e 



d{-An) dAC d ^ dAC' 



de 



dr de 



dr 



'^^dAC\d^_dl i/^dAc: 

^ dr J de de dt \ dr ^ 



(51) 



2 

na 



d{-An) dAC d ^ dAC 



da 



dr da 



dr 



^^^dAC\ds_dfl(^^dAC: 
y dr J da da dt \ dr ^ 
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Similarly, the gauge-invariant Delaunay-type system can be written down as: 

dL _ d ( - AH) _ dAL d / aA£ \ / d A£ \ dg _ dr_ d_ f d AC \ 
H ~ dMo Wm^o\ ^ dr )^\ ^ dr ) dWo ~ SM^ Jt [ ^ dr ) 




dMo ^_ dj-AH) dAC d_ dAC\ ^ dAC.\ dg dr ^ dA£ \ 

dt dL dr dL\ dr J [ dr J dL dL dt\ dr J 

dG _ d{-An) _ dAC d_ dAC \ _ d AC \ dg 

dt duj dr du \ dr J \ dr J duj 

duj _ dj-AH) dAC d / dAC \ 
H ~ dG ^ dr dCy dr J ^ 

dl£ _ d ( - AH) _ dAC d_ f d AC \ / d AC \ dg _ d£ d^ ( dAC \ 
m ~ dn Wdn\ ^ dr ) ^ \ ^ dr ) dQ ~ dnJt\ ^ dr ) 

dn _ _ d ( - AH) dAC d / dAC \ / dAC \ dg 9 AC ' 

Itt ~ dH ^ dr m\ ^ dr ) ^ \ ^ dr J dH ^ dHJt[ ^ dr 

where /i stands for the reduced mass, while 

L = /xV2aV2 ^ G = ii'/'a'/' (l - e'Y^' , H = a'/' (l - e')'^' cosz . (58) 

The symbols /, g now denote the functional dependencies of the gauge, position, and 
velocity upon the Delaunay, not Keplerian elements, and therefore these are functions different 
from /, g used in (46 - 51) where they stood for the dependencies upon the Kepler ele- 
ments. (In Efroimsky (2002a,b) the dependencies /, g upon the Delaunay variables were 
equipped with tilde, to distinguish them from the dependencies upon the Kepler coordinates.) 

To employ the gauge-invariant equations in analytical calculations is a delicate task: one 
should always keep in mind that, in case $ is chosen to depend not only upon time but also 
upon the "constants" (but not upon their derivatives) , the right-hand sides of these equations 
will implicitly contain the first derivatives dCi/dt , and one will have to move them to the left- 
hand sides (like in the transition from (22) to (23)). The choices $ = and $ = —dAC/dr 
are exceptions. (The most general exceptional gauge reads as $ = — dAC/dr -\- rj{t) , 
where r]{t) is an arbitrary function of time.) 

As was expected from (30), both the Lagrange and Delaunay systems simplify in the gauge 
(31). Since for orbital motions we have dH/dp — — dAC/dr , then (31) coincides with (44). 
Hence, the Hamiltonian analysis (34 - 44) explains why it is exactly in the gauge (31) that the 
Delaunay system becomes symplectic. In physicists' parlance, the canonicity condition breaks 
the gauge symmetry by stiffly fixing the gauge (44), gauge that is equivalent, in the orbital 
case, to (31) - phenomenon called "gauge stiffness." The phenomenon may be looked upon 
also from a different angle. Above we emphasised that the gauge freedom imphes essential 
arbitrariness in our choice of the functional form of $(t ; Ci , . . . , Cq) , provided the choice 
does not come into a contradiction with the equations of motion - an important clause that 
shows its relevance in (34 - 44) and (51 - 56): we see that, for example, the Lagrange choice 
^ — (as well as any other choice different from (31)) is incompatible with the canonical 
structure of the equations of motion for the elements. 
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3 A practical example on gauges: 

a satellite orbiting a precessing oblate planet. 

Above we presented the Lagrange- and Delaunay-type planetary equations in the gauge- 
invariant form (i.e., for an arbitrary choice of the gauge function Ci , Cq)) and 
for a generic perturbation AC that may depend not only upon positions but also upon veloc- 
ities and the time. We saw that the disturbing function is the negative Hamiltonian variation 
(which differs from the Lagrangian variation when the perturbation depends on velocities). 
Below, we shall also see that the functional dependence of A7i upon the orbital elements is 
gauge-dependent . 

3.1 The Gauge Freedom and the Freedom of Frame Choice 

In the most compressed form, implementation of the variation-of-constants method in orbital 
mechanics looks like this. A generic solution to the two-body-problem is expressed with 

r = f{C,t) , (59) 




(61) 



c 



- IL J. 

P f 

and is used as an ansatz to describe the perturbed motion: 

r = f{C{t), t) , (62) 

df df dCi ^ 

^ dt dCi dt dt P f dCi dt dt ' ^ ^ 

As can be seen from (63), our choice of a particular gauge is equivalent to a particular way of 
decomposition of the physical motion into a movement with velocity g along the instantaneous 
conic, and a movement caused by the conic's deformation at the rate $ . Beside the fact that 
we decouple the physical velocity r in a certain proportion between these two movements, 
g and it also matters what physical velocity (i.e., velocity relative to what frame) is 
decoupled in this proportion. Thus, the choice of gauge does not exhaust all freedom: one can 
still choose in what frame to write ansatz (62) - one can write it in inertial axes or in some 
accelerated or/and rotating ones. For example, in the case of a satellite orbiting a precessing 
oblate primary it is most convenient to write the ansatz for the planet-related position vector. 

The kinematic formulae (62 - 64) do not yet contain information about our choice of the 
reference system wherein to employ variation of constants. This information shows up only 
when (62) and (64) get inserted into the equation of motion r + {jir/r^) — AF to render 

ag rfQ d# ^ dAC _ d_ ( dAC \ 

dCi dt dt Or dt [ dr ' ^ ^ 



17 



Information about the reference frame, where we employ the method and define the elements 
Ci , is contained in the expression for the perturbing force A/. If the operation is carried 
out in an inertial system, A/ contains only physical forces. If we work in a frame moving 
with a linear acceleration a , then A/ also contains the inertial force — a . In case this 
coordinate frame also rotates relative to inertial ones at a rate /x , then A/ also includes the 
inertial contributions — 2/xxr — /txr — /xx(/ixr). When studying orbits about an 
oblate precessing planet, it is most convenient (though not obligatory) to apply the variation- 
of-parameters method in axes coprecessing with the planet's equator of date: it is in this 
coordinate system that one should write ansatz (62) and decompose r into g and $ . This 
convenient choice of coordinate system will still leave one with the freedom of gauge nomination: 
in the said coordinate system, one will still have to decide what function # to insert in (63). 



3.2 The disturbing function in a frame 

co-precessing with the equator of date 

The equation of motion in the inertial frame is 

where U is the total gravitational potential, and time derivatives in the inertial axes are denoted 
by primes. In a coordinate system precessing at angular rate iJ,{t) , equation (66) becomes: 

r = — z/j, X r — fi X r — fj, X {/J, X r) 

or 



^ ~~dr dr ^tJ'Xr-nxr-nx{nxr) , (67) 

dots standing for time derivatives in the co-precessing frame, and fi being the coordinate 
system's angular velocity relative to an inertial frame. Formula (125) in the Appendix gives 
the expression for /u, in terms of the longitude of the node and the inclination of the equator of 
date relative to that of epoch. The physical (i.e., not associated with inertial forces) potential 
U{r) consists of the (reduced) two-body part Uo{r) = —GM r/r^ and a term AC/(r) 
caused by the planet's oblateness (or, generally, by its triaxiality) . 

To implement variation of the orbital elements defined in this frame, we note that the 
disturbing force on the right-hand side of (67) is generated, according to (65), by 

A£ (r, r,t) ^ - AU{r) + r-{iJ, x r) + ^ (/* x r)-{iJ, x r) . (68) 



Since 



then 



-Q^ = X r , (69) 



•9A£ 

p = r —TTT- = r + fj, X r (70) 
or 
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and, therefore, the corresponding Hamiltonian perturbation reads: 



A7i 



AU + p- (nxr) 



AC + 



1 (dACV 

2 [ dr J 



(71) 



[- AU + {r X p) ■ 1^] = AU - J ■ 1^ , 



with vector 3 = r x p being the sateUite's orbital angular momentum in the inertial frame. 
According to (63) and (70), the momentum can be written as 

p = g + * + /xx/ , (72) 

whence the Hamiltonian perturbation becomes 



AH 



AC + - 



1 fdAcV 



[-AU+ (/ X g) • /z + (# + /z X /) • (/X X /) 



2\ dr J 

This is what one is supposed to plug in (30) or, the same, in (46 - 57). 

3.3 Planetary equations in a precessing frame, 
written in terms of contact elements 

In the preceding subsection we fixed our choice of the frame wherein to describe the orbit. 
By writing the Lagrangian and Hamiltonian variations as (68) and (73), we stated that our 
elements would be defined in the frame coprecessing with the equator. The frame being fixed, 
we are still left with the freedom of gauge choice. As evident from (33) or (46 - 57), the special 
gauge (31) ideally simphfies the planetary equations. Indeed, (31) and (69) together yield 

# = - ^ - fxxr = - nx f , (74) 

wherefrom the Hamiltonian (73) becomes 

ATiM = - [ - AU{f) + M • (/ X g) ] , (75) 
while the planetary equations (30) get the shape 

da d ( - A7i('^°'^*) ) 



or, the same, 

[C.^]^ = ^ [-A[/(/) + M-(/xg)] , (77) 

where / and g stand for the undisturbed (two-body) functional expressions (59) and (60) 
of the position and velocity via the time and the chosen set of orbital elements. Planetary 
equations (76) were obtained with aid of (74), and therefore they render non-osculating orbital 
elements that are called contact elements. This is why we equipped the Hamiltonian (75) with 
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superscript '^(cont)." In distinction from the osculating elements, the contact ones osculate 
in phase space: (72) and (74) entail that p = g . As already mentioned in the end of 
section 1, existence of such elements was pointed out by Goldreich (1965) and Brumbcrg et 
al (1971) long before the concept of gauge freedom was introduced in celestial mechanics. 
Brumberg et al (1971) simply defined these elements by the condition that their insertion 
in g(t ; Ci , . . . , Cq) returns not the perturbed velocity, but the perturbed momentum. 
Goldreich (1965) defined these elements (without calling them "contact") differently. Having 
in mind inertial forces (67), he wrote down the corresponding Hamiltonian (71) and added 
its negative to the disturbing function of the standard planetary equations (without enriching 
the equations with any other terms). Then he noticed that those equations furnished non- 
osculating elements. Now we can easily see that both Goldreich's and Brumberg's definitions 
follow from the gauge choice (31). 

When one chooses the Keplerian parameterisation, then (77) becomes: 



da _ ^ a(-A7i(^""*)) 
dt na dMo 



(78) 



de _l-e^ d{- A7^(-"*)) (i _ d (- AH^ 



cont) 



dt na"^ e dMo na^ e du 



(79) 



du^ ^ -cosz a(-A7i(^°"^)) (i-e^)V^ ^(-AH^-^)) 

dt na?{l — e^yl'^ sinz di na?e de 



di, _ cosi g(_A7i(^°^^)) 1 a(-A7i(^°"*)) 

dt na^(l — e^)^/^ sin i dcu na^(l — e^)^/^ sini 80, 



dn _ 1 a(-A7i(^""*)) 

dt ~ na^l - e2)V2 gin i di ' ^^^^ 



dMo _ _ 1 - a(-A7i(^"*)) _ j_ a(-A7i(<^°^^)) 
dt na^ e de na da 

The above equations implement an interesting pitfall. When describing orbital motion relative 
to a frame coprecessing with the equator of date, it is tempting to derive the Hamiltonian 
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variation caused by the inertial forces, and to simply plug it, with a negative sign, into the dis- 
turbing function. This would entail equations (76 - 83) which, as demonstrated above, belong 
to the non-Lagrangc gauge (31). The elements furnished by these equations are nonosculating, 
so that the conies parameterised by these elements are not tangent to the perturbed trajectory. 
For example, i gives the inclination of the instantaneous non-tangent conic, but differs from 
the real, physical physical (i.e., osculating), inclination of the orbit. This approach - when an 
inertial term is simply added to the disturbing function - was employed by Goldreich (1965), 
Brumberg et al (1971), and Kinoshita (1993), and many others. Goldreich and Brumberg 
noticed that this destroyed osculation, Kinoshita failed to. 

Goldreich (1965) studied how the equinoctial precession of Mars influences the long-term 
evolution of Phobos' and Deimos' orbit inclinations. Goldreich assumed that the elements 
a and e stay constant; he also substituted the Hamiltonian variation (75) with its orbital 
average, which made his planetary equations render the secular parts of the elements. He 
assumed that the averaged physical term ( AU ) is only due to the primary's oblateness: 

, ^ -rr^ J2 93 cos^i — 1 

p being the mean equatorial radius of the planet,^ and n being the satellite's mean motion. 
To simplify the inertial term, Goldreich employed the well known formula 



r X = ^^Gma (1 - e"^) w , (85) 

where 

w = Xi sin i sin Q — X2 sin i cos -|- X3 cos i (86) 

is a unit vector normal to the instantaneous ellipse, expressed through unit vectors Xi, X2) ^3 
associated with the co-precessing frame xj, X2, X3 (axes Xi and X2 lying in the planet's 
equatorial plane of date, and Xi pointing along the fiducial line wherefrom the longitude of 
the ascending node of the satellite orbit, fl , is measured). This resulted in 

(A7i(-"*)) = - [-{AU) + (/i-(/xg))] = 

GmJ2 pI 3 cos^i - 1 /— — . . . „ . . „ , 

: — — — \ G ma [1 — e^) ( pi sim smiZ — p2 sim cossZ + ps cosi) , 

all letters now standing not for the appropriate variables but for their orbital averages. Sub- 
stitution of this averaged Hamiltonian in (81 - 82) lead Goldreich, in assumption that both 
|/Lt|/ ( J2 sin i ) and |ju|/ ( n J2 sin i ) are much less than unity, to the following system: 

dQ 3 , /Pe\^ cosi 

— nJ2 ( — ] o , (88 

dt 2 ' \ aj (1- e^? ^ ' 



di 

— ^ — Pi cos ft — p2 sinfi , (89) 

CLv 

^ Goldreich used the nonsphericity parameter J = (3/2) {p^/ pf J2 , where Pe is the mean equatorial radius. 
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whose solution, 



i = - — cos[ - X {t - to) + ^o] + — sm[ - X {t - to) + ^o] + io , 
X X 



(90) 



^ , \ ^ 1 3^/Pe\^ COSZ 

it ^ -X [t - to) + ^lo where x = 2'^ ) 7]~Z 



n\2 ' 



tells us that in the course of equinoctial precession the satellite inclination oscillates about io ■ 
Goldreich (1965) noticed that his i and the other elements were not osculating, but he 
assumed that their secular parts would differ from those of the osculating ones only in the 
orders higher than 0(|/x|) . Below we shall probe the applicability limits for this assumption. 
(See the end of subsection 3.5.) 



3.4 Planetary equations in a precessing frame, 
in terms of osculating elements 

When one introduces elements in the precessing frame and also demands that they osculate in 
this frame (i.e., obey the Lagrange constraint $ = 0), the Hamiltonian variation reads:^*^ 

A7^(--) = -[-At/ + ^.(/xg) + (mx/)-(/xx/)] , (91) 

while equation (30) becomes: 

dCi dAH^"'^^ 



[Cji Ci 



dt dCr, 



(92) 



To ease the comparison of this equation with (77), it is convenient to split the expression (91) 
for ATiS°^'^'> into two parts: 

^^(cont) ^ _ [ R^^^^^^^f^ t) + tX-{fXg)] (93) 

and 

-(/xx/)-(yLix/) , (94) 
and then to group the latter part with the last term on the right-hand side of (35): 

da d AH^'^*^ 



[Cn Ci 



dt da 



n 



(96) 



10 Both ah'™"*) and AW(°'"=) arc equal to - [ - A[/(/, t) + /x • J] = - [ - AC/(/, t) + /x • (/ x p)] . 
However, the canonical momentum now is different from g and reads as: P = g + (/^ x /) . Hence, the 
functional forms of AH^°^''\f , p) and AH^''°'"'\f , p) are different, though their values coincide. 
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Comparison of this analytical theory with a straightforward numerical integration^^ has con- 
firmed that the 0(|/xp) term in (95) may be neglected over time scales of, at least, hundreds 
of millions of years. In this approximation there is no the difference between ATi^™'**'' and 
^-^(osc) ^ shall write down the equations as: 



[Cji Ci 



da 

dt 



+ M 



X g - / X 



dCn, 



- A- / X 



dCn, 



(96) 



For Ci chosen as the Kepler elements, inversion of the Lagrange brackets in (90) will yield 
the following Lagrange-type system: 



da 
'di 



2 

na 



cont) 



dMo 



A- / X 



or 

dMo, 



(97) 



de 
'di 



na^ e 



d 



dMo 



/ X 



Of ' 
dMo, 



(98) 



(1 - e') 
na?' e 



2\l/2 



dijj 



^ (df . 
+ — xg-/x — 
\duj duj , 



dw — cos i 

dt na^(l — e^)^/^ sin i 



+ 



M _g2Nl/2 



na^ e 



d{-^H^ (df ^ dg\ . / dp 

de \ de de \ de . 



di 



cos I 



dt na^ (1 — e^y/'^ sin i 



9(-A7f(-^)) (df ag\ . / dr 

du) \duj duj \ du) , 



na^ (1 — e^y/'^ sin i 



Credit for this comparison goes to Pini Gurfil and Valery Lainey. 
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dn 

'dt 



1 



na^ (1 - e2)V2 



sm I 



d ( - AH^'^""*) ) 



di 



dMp 
If 



1 - 

na^ e 



a f - A7^(^°"*)) (df ds' 

^ + ^.M-xg-/x/ 

oe \ oe oe , 



2 

na 



5a 



(102) 



terms /x • ( (df/dMo) x g — (dg/dMo) x / ) being omitted in (97 - 98), because these terms 
vanish identically (see the Appendix to Efroimsky (2004) ). 



3.5 Comparison of calculations performed in the two above gauges 

Simply from looking at (76 - 83) and (96 - 102) we notice that the difference in orbit descriptions 
performed in the two gauges emerges already in the first order of the precession rate /x and 
in the first order of fi . 

Calculation of the /it- and /t-dependent terms emerging in (97 - 102) takes more than 20 
pages of algebra. The resulting expressions are published in Efroimsky (2005a), their detailed 
derivation being available in web-archive preprint Efroimsky (2004). As an illustration, we 
present a couple of expressions: 



/ X ^ 

di 



2n2 



(1 - e^) 



(1 + e cosily 



{ fii [— cosfl sm{u! + u) — sin fl cos{u! + u) cos i ] sm{u! + u) 



+ fi2 [— sinfl sin(a; + u) + cos Q cos(a; + u) cos i ] sin(a; + u) 



+ fis sin(a; + cos(a; + u) sin i } 



(103) 



/ 9f 



na?' (3 e + 2 cos v -\- cos v) 
(1 + e cos I/) \/l — 



(104) 



24 



V denoting the true anomaly. The fact that almost none of these terms vanish reveals that 
equations (76 - 83) and (96 - 102) may yield very different results, i.e., that the contact elements 
may differ from their osculating counterparts already in the first order of /x . 

Luckily, in the practical situations we need not the elements per se but their secular parts. 
To calculate these, one can substitute both the Hamiltonian variation and the and /it- 
dependent terms with their orbital averages^^ calculated through 

(•••) ^ ^ ,J / ■■■ 7^-^ ^- (105) 

2 TT Jo (1 + e cost/) 

The situation might simplify very considerably if we could also assume that the precession rate 
fi stays constant. Then in equations (97 - 102), we would assume /x = const and proceed 
with averaging the expressions ( {df /dCj) x g — / x (dg/dCj) ) only (while all the terms 
with jj, will now vanish). 

Averaging of the said terms is lengthy and is presented in the Appendix to Efroimsky 
(2004). All in all, we get, for constant /i : 



M-((^xg-/x^))=/z-(:^xg-/x^)=^/.^ V ^^^^^^^ — ' (^°^) 





^ ■ ^ ^ ~ ^) ^ ^ ° ' Cj = e,n,uj,t,M, 

Since the orbital averages (107) vanish, then e will, along with a , stay constant for as long 
as our approximation remains valid. Besides, no trace of /j, will be left in the equations for 

and i . This means that, in the assumed approximation and under the extra assumption 
of constant fi , the afore quoted analysis (84 - 90), offered by Goldreich (1965), will remain 
valid at time scales which are not too long. 

In the realistic case of time-dependent precession, the averages of terms containing n 
and /i do not vanish (except for /x- ( (df/dMo) x g — / x (dg/dMo) ) , which is identically 
nil). These terms show up in all equations (except in that for a ) and influence the motion. 
Integration that includes these terms gives results very close to the Goldreich approximation 
(approximation (90) that neglects the said terms and approximates the secular parts of the 
nonosculating elements with those of their osculating counterparts). However, this agreement 
takes place only at time scales of order millions to dozens of millions of years. At larger time 
scales, differences begin to accumulate (Lainey et al 2005). 

In real life, the equinoctial-precession rate of the planet, fi , is not constant. Since the 
equinoctial precession is caused by the solar torque acting on the oblate planet, this precession 
is regulated by the relative location and orientation of the Sun and the planetary equator. This 
is why /Lt of a planet depends upon this planet's orbit precession caused by the pull from the 
other planets. This dependence is described by a simple model developed by Colombo (1966). 

Mathematically, this procedure is, to say the least, not rigorous. In practical calculations it works well, at 
least over not too long time scales. 



(107) 
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4 Conclusions: how we benefit from the gauge freedom. 



In this article we gave a review of the gauge concept in orbital and attitude dynamics. Essen- 
tially, this is the freedom of choosing nonosculating orbital (or rotational) elements, i.e., the 
freedom of making them deviate from osculation in a known, prescribed, manner. 

The advantage of elements introduced in a nontrivial gauge is that in certain 
situations the choice of such elements considerably simplifies the mathematical 
description of orbital and attitude problems. One example of such simphfication is the 
Goldreich (1965) approximation (90) for sateUite orbiting a precessing oblate planet. Although 
performed in terms of non-osculating elements, Goldreich's calculation has the advantage of 
mathematical simplicity. Most importantly, later studies (Efroimsky 2005a,b) have confirmed 
that Goldreich's results, obtained for nonosculating elements, serves as a very good approxi- 
mation for the osculating elements. To be more exact, the secular parts of these nonosculating 
elements coincide, in the first order over the precession-caused perturbation, with those of their 
osculating counterparts, the difference accumulating only at very long time scales - see the end 
of section 3 above. A comprehensive investigation into this topic, with the relevant numerics, 
will be presented in Lainey et al (2005). 

On the other hand, neglect of the gauge freedom may sometimes produce camou- 
flaged pitfalls caused by the fact that nonosculating elements lack evident physical 
meaning. For example, the nonosculating "inclination" does not coincide with the real, physi- 
cal inclination of the orbit. This happens because nonosculating elements parameterise instan- 
taneous conies nontangent to the orbit. Similarly, nonosculating Andoyer elements L , G , H 
are no longer the same projections of the angular momentum as their osculating counterparts. 



Appendix 1. Mathematical formalities: 

Orbital dynamics in the normal form of Cauchy 

Let us cast the perturbed equation 

r ^ F + Af ^ - + Af (108) 

into the normal form of Cauchy: 

r = V , (109) 
V = - ^"^ + Af{r{t,C,,...,C,), v{t,C,,...,C,),t) . (110) 

Insertion of our ansatz 

r = fit, Ci(t), Ce{t)) , (111) 
will make (109) equivalent to 

V = ^ + E 9^ . (112) 
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The function / is, by definition, the generic solution to the unperturbed equation 



r^r ■ ^^^^^ 



This circumstance, along with (112), will transform (109) into 
-f dC, 



^ || Q + * = AF( fit, Ci,...,C6) , g(i, Ci,...,C6) + * ) (114) 



where 



* - E ^ (115) 

is an identity, /(t , Ci , ... , Cq) and g(t , Ci , ... , Cq) = df/dt being known functions. Now 
(114 - 115) make an incomplete system of six first-order equations for nine variables (Ci , ... , Ce , 
$1 , ... , $3) . So one has to impose three arbitrary conditions on C , $ , for example as 

$ = *(t, Ci,...,C6) . (116) 

This will result in a closed system of six equations for six variables Cj : 

5: ^ a = AF( fit, Ci,...,C6) , git, C,,...,C,) + * ) - $ (117) 

$ = $( t , Ci , ... , Ce ) now being some fixed function (gauge). A trivial choice is 
$( t , Ci , ... , Cq ) = , and this is what is normally taken by default. This choice is only 
one out of infinitely many, and often is not optimal. Under an arbitrary, nonzero, choice of the 
function $( t , Ci , ... , ) , the system (117 - 118) will have some different solution Cj(t) . 
To get the appropriate solution for the Cartesian components of the position and velocity, one 
will have to use formulae 

r - fit, Ci,...,Ce) , (119) 
f = V = g(i, Ci, Ce) + *(t, Ci, Ce) , (120) 



Generally, $ may depend also upon the variables' time derivatives of all orders: ^{t; Ci , Ci , Ci , ...). 
This will give birth to higher time derivatives of C in subsequent developments and will require additional 
initial conditions, beyond those on r and r , to be fixed to close the system. So it is practical to accept (116). 
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Appendix 2. Precession of the equator of date 



relative to the equator of epoch 



The afore introduced vector fx is the precession rate of the equator of date relative to 
the equator of epoch. Let the inertial axes ( X , Y , Z) and the corresponding unit vectors 
( X , Y , Z ) be fixed in space so that X and Y belong to the equator of epoch. A rotation 
within the equator-of-epoch plane by longitude hp , from axis X , will define the line of nodes, 
X . A rotation about this line by an inclination angle Ip will give us the planetary equator 
of date. The line of nodes x , along with axis y naturally chosen within the equator-of-date 
plane, and with axis z orthogonal to this plane, will constitute the precessing coordinate 
system, with the appropriate basis denoted by ( x , y , z ) . 

In the inertial basis ( X , Y , Z ) , the direction to the North Pole of date is given by 



z = 



T 

sin Ip sin hp , — sin Ip cos hp , cos Ip ) 



while the total angular velocity reads: 



(inertial) 
total 



(inertial) 



(121) 



(122) 



the first term denoting the rotation about the precessing axis z , and the second term being 
the precession rate of z relative to the inertial frame ( X , Y , Z ) . This precession rate is 
given by 



(inertial) 



— (^ip cos hp , ip sin hp , hp^ 



(123) 



because this expression satisfies x z = z . 

In a frame co-precessing with the equator of date, the precession rate will be represented 
by vector 



(inertial) 



(124) 



where the matrix of rotation from the equator of epoch to that of date (i.e., from the inertial 
frame to the precessing one) is given by 



cos hr, 



COS Ip sin hp 



sin Ip sin hp 



sin hr, 



cos/p sin hp 



sin Ip sin hp 



sin In 



COS/r, 



From here we get the components of the precession rate, as seen in the co-precessing coordinate 
frame (a; , y , z) : 



/X = ( //I , //2 , A«3 ) = ( 4 ' hp sin Ip , hp cos Ip ) 



(125) 
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